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In  this  thesis,  I have  analyzed  a notch  in  a leuninated 
composite  plate  by  nvunerically  modelling  the  growth  of  the 
crack  tip  damage  zone.  The  pseudo-crack  represented 
by  the  damage  growth  is  correlated  with  the  actual  lamina  sub- 
crack growth  which  is  shown  in  photographs.  The  growth  of 
the  damage  zone  is  shown  through  a series  of  plots.  Using 
three  different  methods,  I was  able  to  predict  the  fracture 
strength.  It  is  hoped  that  the  results  presented  will  increase 
the.  understanding  of  fracture  in  laminated,  composite  structures. 
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A NUMERICAL  ANALYSIS 


1 


OF  FRACTURE  IN  A LAMINATED, 

FIBROUS  COMPOSITE  PLATE 

I.  Introduction 

This  thesis  is  specifically  concerned  with  the  use  of 
numerical  methods  to  model  a crack  in  a laminated  composite 
plate.  The  crack  is  simulated  by  a through-the-thickness,  fin- 
ite width,  center  notch  oriented  normal  to  the  loading  axis. 
Conventional  finite  element  analysis  and  classical  laminated 
plate  theory  are  used  in  the  numerical  model.  These  techniques 
are  coupled  with  the  use  of  composite  strength  theory  and  in- 
cremental loading  to  follow  the  growth  and  development  of  the 
dcunage  zone  at  the  crack  tip.  There  are  a variety  of  reasons 
why  modelling  should  be  addressed,  but  the  primary  motivating 
factors  are  money  and  safety. 

Background 

Recent  fuel  price  increases  have  been  the  motivation  for 
increased  research  into  fuel-conservation  technology.  Signi- 
f,icant  fuel  savings  can  be  realized  by  lighter  weight  vehicles, 
and  one  of  the  most  promising  ways  of  reducing  the  structural 
weight  of  aerospace  vehicles  is  to  use  high  strength  composites 
in  primary  structures  [1] . As  recognized  by  current  Air  Force 


policy,  fracture  considerations  are  important  in  aircraft 
design  and  are  especially  important  for  the  primary  structures 
of  aircraft  [2].  Although  current  Air  Force  fracture  design 
requirements  only  apply  to  metal  aircraft  structures,  as  com- 
posites become  more  widely  applied,  fracture  mechanics  con- 
siderations in  the  design  of  composite  components  will  surely 
become  mamdatory  in  the  interest  of  safety. 

Linear  elastic  fracture  mechanics  (LEFM)  theory  has  been 
developed  for  describing  the  behavior  of  brittle,  homogeneous, 
isotropic  materials.  LEFM  can  be  used  to  describe  the  three  (3) 
basic  modes  of  crack  extension,  shown  in  Figure  1 [3].  The  crack  follows 
three  stages  of  growth.  First  a crack  initiates  from  an  exist- 
ing flaw;  it  then  propagates  in  a stcUDle,  usually  slow  manner, 
amd  last  comes  unstable  extension  and  structural  failure  [4] . 

The  crack  will  propagate  in  the  direction  along  which  the  elas- 
tic energy  release  rate  per  unit  crack  extension  will  be  maximum. 
In  Mode  I extension,  this  direction  is  perpendicular  to  the  di- 
rection of  greatest  local  tension  [5] . 

In  LEFM,  crack  propagation  is  explained  by  an  energy  balance 


extension,  and  the  material  in  which  the  crack  is  propagating 


[51 . The  basic  theory  assumes  that  the  entire  body  containing 


the  crack  is  elastic,  when  in  actual  materials  the  stress  sin- 


gularity at  the  crack  tip  produces  plastic  flow. 


Figure  1.  Modes  of  Crack  Extension 


LEPM  can  be  modified  to  account  for  plasticity  if  the 


plastic  region  is  small  compared  to  the  crack  length  [7] . The 


plastic  deformation  at  the  crack  tip  effectively  blunts  the  tip 


and  hence  makes  the  material  tougher  [8] . Several  models  are 


available  to  accoimt  for  the  effect  of  plasticity,  but  the  sim- 
plest method  is  to  assume  a plane  stress  solution  and  then 


modify  the  crack  length  by  a parameter  determined  by  the  region 


where  yielding  has  occurred  [8] . 


If  the  small  area  at  the  crack  tip  in  fibrous  composites 


could  be  accurately  modelled  as  homogeneous  and  anisotropic. 


many  of  the  relations  from  LEFM  could  be  diiectly  applied  since 
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none  of  the  basic  principles  used  in  fracture  mechanics  would 
be  be  violated  [9] . But,  it  has  been  shown  that  fracture  is 
very  sensitive  to  the  local  properties  at  the  crack  tip  [6] ; 
therefore,  the  presence  of  two  phases  complicates  the 
fracture  process  in  composites.  Insight  into  composite  frac- 
ture can  be  gained  by  examining  the  observed  phenomena  of 
composite  fracture. 

Most  of  the  composite  fracture  experiments  have  emphasized 
Mode  I loading  [10] . The  cracks  have  either  been  parallel  or 
perpendicular  to  fibers  in  unidirectional  composites,  or  align- 
ed with  a material  axis  in  composite  laminates  [6] . The  ob- 
servations made  in  these  experiments  can  best  be  understood  by 
relating  them  to  the  three  stages  of  crack  growth. 

Microcracks  initiate  in  the  matrix  almost  from  the  onset 
of  loading.  The  microcrack  initiation  sites  are  flaws  which 
are  introduced  during  the  fabrication  process  [11]  or  are  in- 
duced by  stress  concentration,  battle  damage,  fatigue,  or  tran- 
sient high  loading  conditions  [12] . 

Due  to  the  microscopic  heterogeneity  at  the  crack  tip, 
cracks  in  composites  do  not  propagate  in  the  same  manner  as  they 
do  in  isotropic  materials  [13] . After  initiation,  there  is 
normally  no  visible  self-similar  crack  growth  [14] . Rather, 
the  crack  propagates  by  developing  a network  of  microcracks 
[10] . Due  to  the  difference  in  material  properties  between 
the  fiber  and  the  matrix,  three  types  of  crack  propagation 
can  occur  when  the  microcracks  reach  the  fiber  matrix  interface. 
The  cracks  can  be  reflected  back  into  the  matrix;  they  can 
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continue  to  travel  directly  through  the  fiber,  or  they  may 
cause  interfacial  debonding  [13] . Normally,  subcracks  are 
formed  which  extend  parallel  to  the  fibers.  It  has  been 
observed  that  these  subcracks  continue  to  grow,  either  along 
the  interface  or  in  the  matrix, and  they  are  influenced  by  what 
occurs  in  the  neighboring  fibers  [5] . 

In  l2uninates,  these  subcracks  extend  in  each  ply  along 
the  fibers  in  a manner  similar  to  the  way  the  crack  tip 
plastic  zone  grows  in  metals  [4,13].  In  fact,  one  group  of 
experiments  has  shown  that  this  damage  zone  can  be  treated  like 
the  plastic  zone  in  metals  to  obtain  fracture  strength  [13] . 

The  deunage  zone  seems  to  function  like  a plastic  zone  in  that 
the  "yielding"  in  the  composite  serves  to  blvmt  the  crack  tip 
emd  relax  the  stresses  [4] . 

Subcracks  and  damage  zone  growth  ultimately  lead  to  failure 
of  the  structure.  Ultimate  failure  may  either  be  characterized 
by  ply  failure,  where  the  fibers  actually  pull  out  from  the  ma- 
trix, or  break, by  delamination,  or  by  a combination  [15]. 

Thus  at  the  present  time,  the  direct  applicability  of  LEFM 
to  launinated  composite  fracture  cannot  be  assumed  [4,  6,  13]. 

The  presence  of  notch  sensitivity  [4]  and  the  fact  that  crack 
growth  parallel  to  the  fibers  in  unidirectional  composites  can 
be  explained  by  a stress  intensity  factor  [16,  17]  indicate  that 
some  portions  of  presently  developed  isotropic  fracture  theory 
can  be  applied.  The  major  divergences  from  isotropic  fracture 
theory  are  the  growth  of  a damage  zone  as  opposed  to  crack  open- 
ing from  the  crack  tip  and  the  dissimilar  behavior  of  each  ply 
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in  a laminate  to  a given  load. 

The  data  which  will  be  used  to  validate  the  analysis  done 

in  this  thesis  was  generated  by  the  Air  Force  Materials  Lediora- 

tory  16,  10].  During  the  experiments,  notched  composite  plates 

of  Thomel  300  graphite  fibers  in  Narmco  5208  epoxy  were  loaded 

to  failure.  Figure  2.  The  specimens  shown  in  Fig.  2,  contained 

either  a 13mm  or  a 20mm  notch  in  (0,  +45)  and  (0,  +45,  90) 

— s — s 

laminates.  The  results  of  these  experiments  were  used  to  show 
that  the  notched  fracture  strength  could  be  predicted  using 
unnotched  failure  strength  and  the  dimension  of  the  deunage  zone 
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Figure  2.  Specimen  Geometry 


As  part  of  these  experiments,  radiographs  were  taken  at 
various  load  levels.  The  radiographic  image  of  the  crack 
tip  deunage  area  was  enhanced  as  described  in  Ref.  18.  From 


these  radiographs  it  is  possible  to  measure  the  length  of  the 
subcracks  in  each  ply  of  the  laminate. 

This  thesis  will  use  a finite  element  model  of  the  experi- 
mental (0,  +45,  90)  l^uninate  with  a 13mm  (.51  in.)  notch  to 

S 

analyze  the  growth  of  the  crack  tip  damage  zone.  The  crack 
tip  stress  and  displacement  fields,  and  the  dimensions  and  shape 
of  the  damage  zone  in  each  ply  are  obtained  as  output.  This 
output  data  is  compared  with  subcrack  dimensions  obtained  from 
the  radiographs,  isotropic  crack  stress  and  displacement  fields, 
and  with  several  models  used  to  predict  fracture  strength. 

The  following  sections  explain  the  analysis.  First,  the 
theory  behind  the  modelling  will  be  explained.  Next,  the 
numerical  analysis  will  be  covered.  Last,  the  results  and 
conclusions  based  on  comparisons  between  the  analytical  data 
and  experimental  data  will  be  presented. 
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II.  Theory 


The  theoretical  foundation  of  this  analysis  can  be  built 
by  answering  three  questions.  What  are  the  implications  of  a 
two  dimensional  analysis  of  the  crack  tip  damage  zone  in  a 
composite  laminate?  What  composite  strength  theory  is  appro- 
priate for  crack  tip  damage  zone  analysis,  and  how  can  this 
theory  be  applied?  How  can  the  modelled  damage  region  be 
used  to  predict  fracture  strength? 


Implications  of  Two  Dimensional  Analysis 


Several  investigations  have  stated  that  the  stress  field 
around  a crack  in  a composite  plate  is  three  dimensional  [9,6,13] 
The  three  dimensionality  is  caused  by  the  interlaminar  stresses 
at  the  free  edge  of  a crack. 

The  effect  of  the  free  edge  on  the  interlaminar  stresses  in 
an  unnotched  laminate  plate  under  uniaxial  stress  has  been 
partially  explained  by  Pagano  and  Pipes.  Notches,  though,  pre- 
sent a different  type  of  free  edge.  This  difference,  coupled 
with  the  crack  tip  singularity,  make  the  solution  of  the  free 
edge  effects  along  a notch  extremely  difficult.  Therefore,  any 
discussion  on  the  applicability  of  a two  dimensional  analysis 
must  be  made  based  on  the  purpose  of  the  analysis. 

This  analysis  is  concerned  with  modelling  and  analyzing  the 
crack  tip  damage  zone  and  how  it  relates  to  the  fracture  strength 
The  applicability  of  the  two  dimensional  analysis  can  then  be 
assessed  by  comparing  dimensions  of  the  d2unage  zone  and  the  cupeas 


affected  by  interlaminar  stresses,  and  then  determining  whether 
interlaminar  stresses  affect  fracture  strength. 

Mandell,  Wang,  and  McGarry  [4]  performed  a finite  element 
analysis  of  a single  edge  notched  composite  plate  using  a three 
dimensional  finite  element  analysis.  It  should  be  noted  that 
the  interlaminar  stresses  for  an  edge  notch  are  probably  worse 
than  those  for  a center  notch  since  the  edge  notch  is  connected 
to  a stress  free  edge.  Therefore,  if  the  effects  of  the  inter- 
laminar stresses  are  within  the  bounds  of  the  d2unage  zone  for 
this  edge  case,  then  interlaminar  stresses  for  a center  notched 
specimen  should  be  even  less  Significant. 

The  results  of  the  Mandell  analysis  indicate  that  the  ef- 
fect of  the  interlaminar  stresses  was  confined  to  a distance 
equal  to  the  laminate  thickness  along  the  crack  flanks  and 
was  less  ahead  of  the  crack  tip.  Also,  excunination  of  the 
isostress  plots  reveals  that  the  area  where  the  interlaminar 
stresses  were  of  the  same  magnitude  as  the  planar  stresses  was 
much  smaller.  Since  previous  studies  [10,  13,  14,  4]  indicated 
that  the  damage  zone  was  confined  to  the  area  immediately  at 
the  crack  tip  or  ahead  of  it,  and  the  damage  zone  at  failure 
was  much  greater  than  the  laminate  thickness,  the  interlaminar 
stresses  should  not  affect  the  boundary  of  the  damage  zone. 

In  reference  [19],  the  effects  of  interleuninar  stresses  on 
the  fracture  strength  of  center  notched  laminated  plates  was 
studied  experimentally.  This  study  showed  that  the  interleiminar 
stresses  hal  no  effect  on  the  fracture  strength  of  a center 
notched  plate. 
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I In  light  of  previous  work,  it  seems  feasible  to  use  a two 

dimensional  analysis  to  model  the  boundary  of  the  damage  zone 
in  a center  notched  Icuninated  plate.  Furthermore,  the  stress 

I ! 

field  value  outside  a region  equal  to  a laminate  thickness  from 
the  crack  should  be  accurate. 


j 


Strength  Determination 

In  general  for  any  material,  strength  is  a measure  of  the 
Ability  to  deform  without  sustaining  irreversible  deunage.  In 
homogeneous  isotropic  materials,  this  ability  is  measured  by 
the  yield  criterion,  where  the  yield  strength  is  the  point 
where  the  material  ceases  to  act  elastically.  The  yield  cri- 
terion defines  a hypersurface  against  which  various  loading 
conditions  can  be  evaluated  [20,21]. 

Several  studies  have  modelled  the  behavior  of  a unidirec- 
tional composite  as  linear  to  failure  [9,  22,  23,  24].  Drawing 
a parallel  with  plasticity  theory,  it  can  then  be  surmised  that 
the  combinations  of  stresses  which  represent  lamina  failure 
can  be  represented  by  a hypersurface  in  stress  space  [17] . 
Leunina  failure  is  defined  as  the  Inability  of  the  lamina  to 
carry  stress  in  the  same  manner  as  it  did  in  its  virgin  state. 

After  putting  laminae  together  to  form  a laminate,  the 
behavior  of  the  Izuninate  is  no  longer  linear  to  failure.  In- 
stead it  can  be  assumed  that  failure  of  the  launinate  occurs 
when  all  of  its  constituent  laminae  have  failed  [9,  17]. 

There  are  two  basic  types  of  criteria  which  have  been  de- 
veloped to  describe  the  plastic  yield  surface,  non-interacting  and 
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interacting  [17] . Li  the  nonrinteracting  criteria,  svxdi  as  naxisun.  stress 
or  maxiiiun  strain,  nultiaxial  stress  does  not  affect  uniaxial  strength  [21] . 
Fbr  the  interacting  criteria,  such  as  the  von  Mises  or  Tresca  criteria,  the 
yield  surface  is  determined  by  the  total  stress  tensor.  Since 
this  analysis  deals  with  a high  strength  graphite  epoxy  and 
examination  of  failure  curves  for  high  strength  graphite-epoxy 
specimens  [25]  reveals  a high  degree  of  Interaction,  the  non- 
interacting criteria  will  not  be  used  in  this  analysis. 

There  are  basically  two  interaction  theories  which  have 
been  developed  for  homogeneous,  anisotropic  materials.  One 
is  the  Tsai-Hill  theory  [20,  22],  in  which  the  failure  surface 
is  defined  by 

F(o  -o  + G(a  -o  )*  + H(a  -a  )*  + 

2 3 3 1 1 2 

2Lt  * + 2Mt  * + 2Nt  * =»  1 (1) 

23  31  12 


where  F,  G,  H,  L,  M,  and  N are  failure  strength  par2uneters. 
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The  other  theory  is  the  Tsai-Wu  tensor  theory,  in  which  the 
failure  surface  is  defined  by 


Vi  * * • 


where  F^^  and  F^^^  are  strength  tensors  of  the  second  and  fourth 
order  respectively  [9] . Although  more  general  than  the  Tsai- 
Hlll  criterion,  the  Tsai-Wu  criterion  is  more  complicated  and 
requires  extensive  testing  to  determine  the  values  of  the 
strength  tensors.  Since  this  work  is  analytical  and  further 
testing  is  required  to  determine  the  strength  tensor,  the 


Tsai-Hill  criterion  will  be  used  since  the  failure  strength 

parameters  can  be  determined  from  the  uniaxial  strengths  for 
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o*  a a a*  a* 

+ -i-  + = 1 (3) 

X*  X*  Y*  S* 

where  X is  axial  strength,  Y is  transverse  strength  and  S is 
shear  strength. 

When  the  Tsai-Hill  criterion  is  used,  extreme  caution  is  re- 
quired for  two  reasons  [221 . First,  this  criterion  is  based 
on  a distortional  energy  approach  where  hydrostatic  pressures 
do  not  cause  yielding  [20] . It  has  been  shown  that  composites 
do  yield  under  hydrostatic  loading  due  to  the  misalignment  of 
the  principle  stress  and  principle  strain  directions  [23] . 

Second,  the  criterion,  as  stated,  does  not  account  for  differ- 
ences between  compressive  and  tensile  strengths.  The  effect 
of  these  deviances  can  be  found  by  examining  the  applicability 
of  the  criterion  to  the  material  being  analyzed. 

This  criterion  is  phenonenalogical  in  nature;  therefore,  the- 
applicability  of  the  criterion  must  be  judged  against  its  abil- 
ity to  predict  failure  for  the  particular  laminate  under  study. 

It  will  be  shown  that  the  Tsai-Hill  criterion  can  be  used  to 

predict  uniaxial  failure  of  the  (0,  +45,  90)  graphite-epoxy 

s 

laminate  used  in  this  analysis;  therefore,  it  will  be  assumed 
that  the  problems  arising  out  of  yielding  under  hydrostatic 
loading  do  not  significantly  affect  strength  predictions. 

The  difference  in  tensile  and  compressive  failure  strengths 
can  be  included  by  using  a simple  procedure  [9] . If  the  stress 
field  includes  compressive  stresses,  the  compressive  failure 
strength  is  used  for  that  component  of  the  stress  field,  otherwise 
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the  tensile  failure  strength  is  used  in  the  criterion. 

Before  this  criterion  can  be  applied  to  a laminate,  a 
determination  must  be  made  as  to  how  to  treat  constituent 
lamina  failure  when  it  occurs  prior  to  total  l2uninate  failure. 
There  are  various  methods  to  predict  laminate  failure,  the 
most  appropriate  is  the  use  of  the  individual  failure  char- 
acteristics of  the  l2uninae  to  predict  laminate  failure  through 
progressive  lamina  failure  [9,  17,  24].  As  a leunina  fails, 
the  equivalent  stiffness  for  the  leunina te  is  changed  by  modi- 
fying the  constituent  lamina  stiffness. 

There  are  basically  three  methods  to  modify  the  leuninate 
stiffness  [17] . There  is  the  total  discount  method  where  the 
failed  ply  is  assigned  zero  stiffness  and  strength.  Second, 
there  is  the  mode  limited  discount  method  where  zero  stiffness 
and  strength  are  assigned  to  the  transverse  and  shear  modes  if 
the  failure  occurs  only  due  to  transverse  or  shear  stresses, 
and  the  strength  and  stiffness  in  all  modes  is  discounted  if 
the  failure  is  due  to  longitudinal  stresses.  In  the  third 
method,  the  failed  lamina  is  assigned  residual  properties. 

Since  the  exact  lamina  failure  mechanisms  and  the  post 
failure  performance  of  the  launina  in  a laminate  are  not  com- 
pletely understood,  this  analysis  uses  a bounding  theory  to 
portray  failure  in  the  crack  tip  damage  zone.  It  is  obvious 
that  one  extreme  bound  on  the  damage  zone  can  be  found  if 
all  lamina  are  treated  as  elastic  until  failure  occurs  (elastic 
method),  and  the  other  extreme  bound  is  determined  by  completely 
discounting  all  stiffness  of  a ply  when  the  stress  field  in 


that  ply  exceeds  the  failure  criterion. 

The  appropriateness  of  two  dimensional  analysis,  the 
Tsai-Hill  failure  criterion,  the  progressive  Icunina  failure 
approach  to  predict  laminate  failure,  and  the  boundedness 
assumption  made  in  the  analysis  can  be  illustrated  by  applying 
these  procedures  to  an  unnotched  specimen  under  tension.  The 
properties  of  the  graphite-epoxy  specimen  used  are  shown  in 
Table  I. 


Table  I 

Leunina  Properties 

Elastic 

Constants 

Ultimate 

Strengths 

t . 

.0436 

in. 

*t 

217.6 

ksi 

®11 

15545 

ksi 

*c 

217.6 

ksi 

®22 

1425 

ksi 

^t 

5.8 

ksi 

^12 

903 

ksi 

Y 

c 

35.7 

ksi 

'^12 

.288 

S 

9.9 

ksi 

The  ultimate  strength  of  the  (0,  +45,  90)  Icuninate  was 

S 

54.4  ksi. 

Using  these  properties  and  the  procedures  mentioned  pre- 
viously, the  ultimate  strength  of  the  leiminate  was  predicted. 
The  loading  curve  is  shown  in  Fig.  3,  and  the  resultant  forces 

and  strains  at  the  knees  are  shown  in  Table  II.  The  upper  line 
< 

represents  the  curve  obtained  if  all  lamina  remain  elastic  un- 
til all  lamina  have  failed.  The  middle  curve  uses  the  mode 


limited  discount  method,  and  the  lower  curve  uses  the  total 
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Figure  3.  Loading  Curve  Showing  Difference 
Between  Discounting  Methods. 


1 Table  II  1 

Loading  Curve 

Data 

Ply  Which 

Elastic  Method 

Mode  Limited 

Complete  Discount 
Method 

(lb. /In.) 

N 

(lb. /In.) 

X ^ 

(lb. /In.) 

o 

o 

1230  4.4x10“^ 

1230  4.4xl0"^ 

1230  4.4xl0"^ 

+45® 

0® 

3900  14.0x10"^ 

1620  6.1x10“^ 

-3 

3070  14.0x10 

1470  6.2x10"^ 

-3 

2370  14.0x10 

discount  mevnod. 

As  can  be  seen,  the  elastic  method  is  an  upper  bound  on 
failure,  and  the  total  discount  method  does  provide  a lower 
bound.  The  closeness  of  the  lower  bound  to  the  actual  strength 
suggests  that  the  ability  of  the  composite  to  handle  loads  is 
severely  degraded  when  transverse  or  shear  failure  occurs  as  is 
assumed  in  the  total  discount  method. 

In  summary,  failure  in  the  damage  zone  will  be  bounded  using 
the  elastic  method  and  the  complete  discount  method.  Lamina 
failure  will  be  defined  using  the  Tsai-Hill  criterion. 

Application  of  Numerical  Analysis  Results 

The  goal  of  a fracture  analysis  is  to  predict  in  what  fashion 
and  when  catastrophic  failure  will  occur.  The  manner  of  failure 
can  be  correlated  with  either  crack  propagation  or  the  growth 
of  a damage  zone  at  the  crack  tip.  The  "when"  of  ultimate  fail- 
ure is  explained  by  the  fracture  stress  if  strain  rates  and  ■" 
other  time  related  phenomena  are  ignored.  The  phenomena  cor- 
related with  these  two  questions  will  be  discussed  in  the  fol- 
lowing paragraphs. 

Since  crack  extension  did  not  occur  for  the  tests  being 
analyzed,  the  investigation  must  concentrate  on  explaining  the 
growth  of  the  damage  zone.  Two  specific  phenomena  are  related 
to  the  growth  of  the  damage  zone.  The  unobservable  phenomenon 
is  the  growth  of  the  failure  area  around  the  crack  tip.  This 
phenomenon  is  described  by  the  dimension,  c,  of  the  failure 
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area  or  damage  zone  where  all  plies  have  failed,  the  laminate 
damage  zone.  In  this  laminate,  this  zone  coincides  with  the 
dimension  of  the  0**  damage  zone  measured  col  inear  with  the 
original  notch  and  can  only  be  found  through  numerical  analysis. 
The  observable  occurrence  is  the  growth  of  subcracks  along 
ply  fibers  [4,  13]. 

It  is  questionable  whether  the  growth  of  subcracks  and 

the  fracture  stress  are  directly  related  [4] . It  has  been 

shown,  though,  that  the  length  of  the  subcracks  is  proportional 

2 

to  the  opening  mode  stress  intensity  factor,  squared,  [4] . 

The  average  length  of  the  ply  subcracks  can  be  obtained  from 

radiographs  of  the  crack  tip  damage  zone,  as  can  be  seen  in 

2 

Fig.  4,  5,  6,  7,  8,  and  9.  Since  is  proportional  to  the 
strain  energy  release  rate,  [5,  26],  it  is  also  possible  to 
relate  to  the  ply  subcrack  length.  Values  of  and  M can 
be  found  using  the  finite  element  analysis  where  the  crack  ex- 
tension is  equal  to  the  growth  of  the  laminate  deunage  zone. 

This  correlation  indirectly  relates  the  growth  of  the  modelled 
damage  and  the  growth  of  the  ply  subcracks,  which  can  only  be 
measured  through  experiments. 

Since  the  damage  zone  in  laminates  is  analogous  to  the 
crack  tip  plastic  zone  in  metals,  the  following  expression  can 
be  used  to  calculate  in  a center  notched,  finite  width  plate 
[27]. 


K = a (w  tan  Tra) 
w 
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Figure  4, 


Enhanced  Radiographic  Image  of  Ply  Subcracks  at 
50%  of  the  Experimental  Fracture  Load,  Approxi- 
mately Five  Times  Actual  Size. 
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Figure  5.  Enhanced  Radiographic  Image  of  Ply  Subcracks  at 
60%  of  the  Experimental  Fracture  Load,  Approxi- 
mately Five  Times  Actual  Size. 
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Figure  6. 


Enhanced  Radiographic  Image  of  Ply  Subcracks  at 
70%  of  the  Experimental  Fracture  Load,  Approxi- 
mately Five  Times  Actual  Size. 


Enhanced  Radiographic  Image  of  Ply  Subcracks  at 
95%  of  the  Experimental  Fracture  Load,  Approxi- 
mately Five  Times  Actual  Size. 


In  this  expression,  a is  the  applied  stress,  w is  the  plate 
width,  and  a is  the  effective  crack  half-length.  The  effec- 
tive crack  half-length  is  the  sum  of  the  original  crack  half- 
length,  a^,  and  the  size  of  the  dcunage  zone,  c,  measured  co- 
linear  with  the  original  notch,  see  Figure  2. 


a = a + c 
o 


With  this  relation  and  the  results  from  the  finite  element 
analysis,  it  is  possible  to  calculate  for  each  loading  in- 
crement. Since  the  subcrack  length  can  be  related  to  the 

2 

loading  increment,  it  is  then  possible  to  correlate  and 
the  ply  subcrack  lengths. 

Various  studies  have  used  a finite  element  analysis  to 
determine  the  strain  energy  release  rate  [6,  28].  It  has 
been  shown  that  [5,  6] 


2/  - 2-  ^ 
2 da 


HP 

where  P is  the  applied  load  and  is  the  rate  of  change  of 

da 

structural  compliance  with  crack  extension.  From  the  finite 


element  analysis  in  which  element  plies  are  removed  as  they 
fail,  load-displacement  diagrams  are  obtained  similar  to 

I 

those  shown  in  references  [14]  and  [27]  and  in  Fig.  10.  As  the 
damage  zone  increases  the  compliance  changes.  The  compliance 
can  be  calculated  by  assuming  the  displacement  at  the  load 
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application  point  returns  to  zero  as  the  load  is  relaxed, 
line  AO.  The  inverse  of  the  slope  of  this  line  is  the  new 
compliance,  Cj*  The  change  in  crack  extension  is  the  dif- 
ference between  the  equivalent  crack  length  2a2  and  the  pre- 
vious equivalent  crack  length  2a The  derivative  can  then 
be  approximated  by 


^2  ~ ^1 
2a2~  2a2^ 


The  equivalent  strain  energy  release  rate  is  approximated  by 


^2  - ^1 
2a2-  2ai^ 


The  second  goal,  of  predicting  when  fracture  will  occur, 
can  be  approached  in  three  manners.  Two  of  these,  the  use  of 
a critical  stress  intensity  factor  and  instability  analysis 
are  classical  in  nature,  and  the  third  is  introduced,  for  the 
first  time,  in  this  thesis.  The  third  method  involves  the 
relation  between  combinations  of  load  and  load  bearing  area 
which  result  in  stresses  above  the  ultimate  stress. 

In  a previous  work  [13],  Hahn  predicted  that  K /a  is 

c o 

w w 

between  .7393  in.  ^ and  .7663  in.  for  the  specimen  being 
analyzed,  where  is  the  critical  opening  mode  stress  in- 
tensity factor  and  is  the  unnotched  tensile  strength. 

Using  these  values  and  the  unnotched  tensile  strength  of 


40.22  ksi-in.**  < K < 41,67  ksi-in.**  (9) 

c 


Using  equations  (4) , (5) , and  (9)  and  the  dimensions  of  the 
damage  zone  obtained  from  the  finite  element  analysis,  the 
failure  load  can  be  predicted. 

The  instability  method  of  predicting  fracture  strength 
is  the  simplest  to  apply.  The  load-displacement  curve  for 
the  sequence  of  loading  where  element  plies  within  the  fail- 
ure region  are  removed  is  plotted.  Fig.  11.  The  increment 
of  load  which  causes  the  curve  to  transition  from  the  non- 
linear region,  II,  to  the  flat  region.  III,  is  defined  as 
the  instability  load,  and  the  preceding  load  is  taken  as  the 
fracture  strength. 


Figure  11.  Instability  Load  Displacement  Diagraun 
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The  last  method  is  based  on  logic  similar  to  that  pre- 
sented in  an  article  by  Nuismer  and  Miitney  [29] . The  average 
stress  criteria  predicts  failure  when  the  average  value  of 


stress,  a , over  some  fixed  distance,  d , ahead  of  the  crack 

a O 


first  reaches  the  unnotched  tensile  strength,  a^.  In  equa- 


tion form,  failure  occurs  when 
1 


a+d 

s ° 


a o^(0,y)dy  = 


(10) 


where  d^  is  the  fixed  distance  ahead  of  the  crack,  a is  the 


crack  length,  and  is  the  perpendicular,  normal  stress 


component  ahead  of  the  crack  tip.  Clearly  the  average  stress 
is  represented  by  the  left  side  of  equation  (10) . For  a center 
notched  plate.  Fig.  2,  equation  (10)  must  hold  at  both  ends  of 
the  crack,  and  the  normal  stress  perpendicular  to  the  crack 


flemks  must  be  zero;  therefore, 

a+d 


(11) 


-a-d. 


The  equation  as  presented  in  the  referenced  paper  predicted 
the  notched  strength  for  an  infinite  plate  with  a center  notch. 


For  the  finite  plate  being  analyzed,  it  is  assumed  that  d^  is 


the  distance  from  the  boundary  of  the  damage  zone  to  the  plate 
edge. 


2d^  = W - 2a 
o 


(12) 


1 


Further,  if  it  is  assumed  that  the  effective  crack  length  is 

the  sum  of  the  initial  notch  length,  a , and  the  dimension  of 

o 
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the  damage  zone  colinear  with  the  original  notch,  c,  such  that 

(13) 


a a a + c 
o 


From  Fig.  2,  it  can  be  seen  that 


- w , 

“5“  * ’^"*^0 

(14) 

w * a+d 

T"  ° 

(15) 

2d^  = w-2a  -2c 
o o 

(16) 

Incorporating  equations  (14),  (15),  and  (16)  into  equation  (11) 

w 


- v-2a'-  5c 


a w- 


dy 


(17) 


w 

-I 


From  classical  laminated  plate  theory 

f 

N = J a dz 
X y.. 


* -t/2  * 


(18) 


where  is  the  resultant  force  per  unit  length,  and  t is  the 


laminate  thickness.  Now  integrating  equation  (17)  across  the 
thickness  of  the  plate  the  following  is  obtained 


,t/2 

J o dz  = 


1 J 

J w-2a  -20  J o 


-t/2 


-t/2 


,dydz  (19) 


assuming  sufficient  continuity  in  and  that  W,  a^,  and  c 
are  constant  through  the  thickness,  equation  (19)  is  rewritten 
as 


rt/2 

J a dz  = 
-t/2  ^ 


w 


.t/2 


= — = s—  f 2 r cr  dzdy 

w-2a^-2c  X 


(20) 


By  the  definition  of  the  averaging  process,  a is  constant 

O 


through  the  thickness,  therefore 
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t/2 

oAz 

a 


(21) 


From  equations  (18)  and  (19) 


to. 


w-2a  -2c 
o 


/ 


w 

2 


-w 

2 


N^dy 


(22) 


From  energy  considerations,  the  integral  of  the  resultant 
force  in  the  x direction  integrated  over  a line  perpencicular 
to  the  x-axis  must  equal  the  applied  load  in  the  x direction. 


^ax' 


Substituting  equation  (23)  into  (22) 

^°a  ~ w-2a^-2c  ^ax 
or 


o 

a 


1 

t (w-2a^-2c) 


P 

ax 


(23) 


(24) 


(25) 


Defining  the  quantity,  t(w-2a^-2c),  as  the  remaining  load 
bearing  area,  and  substituting  equation  (25)  and  into 

equation  (10) , the  following  condition  will  define  fracture 
strength 


(26) 


To  find  the  fracture  strength  for  a particular  material, 

equation  (26)  is  applied  to  a load  versus  load-bearing  area 
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diagram,  Fig.  12.  The  boundary  of  the  failure  region  is  the 


line  defined  by 


As  damage  progresses  from  the  crack  tip,  the  load  bearing  area 
decreases.  If  load  and  remaining  load  bearing  area  are  plotted 

as  shown  in  Fig. 12,  the  predicted  fracture  strength  is  the  in- 

* 

tercept  between  this  curve  and  the  boundary  of  the  failure  re- 
gion. Referring  to  the  boundedness  argument  presented  in  the 
Strerwth  Determination , the  elastic  method  yields  an  upper 


bound  on  strength,  and  the  complete  discount  method  yields  a 
lower  bound 


Actual 
Fracture 
Load  ^ 


\ \ 

Progressive 

Fciilure 

Model 

A (in.^) 


^Elastic 

Model 


Figure  12.  Applied  Load  Versus  Load  Bearing  Area. 

This  is  the  last  section  of  the  theory  chapter.  Now  that 
the  theory  behind  the  application  of  the  numerical  results  has 
been  built,  it  is  necessary  to  establish  the  foundation  of  the 
finite  element  modeling  which  was  done  to  obtain  the  results. 
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III.  Numerical  Analysis  Description 

The  analysis  generated  in  this  thesis  is  based  on  data 
obtained  through  a finite  element  method.  In  this  section, 
the  finite  element  method  used  and  the  method  employed  to  de- 
termine the  size  of  the  d^unage  zone  in  each  ply  will  be  dis- 
cussed. 

Finite  Element  Model 

To  apply  some  of  the  methods  of  LEFM  or  to  obtain  the  exact 
stresses  and  displacements  around  a crack  tip,  an  exact  elas- 
ticity solution  is  required  [30] . Since  exact  solutions  for 
composite  problems  are  difficult,  if  not  impossible,  to  obtain, 
it  is  appropriate  to  use  approximate  numerical  methods  to  an- 
alyze composite  fracture. 

The  finite  element  progreun  used  in  this  thesis  was  developed 
at  the  Air  Force  Flight  Dynamics  Laboratory  [31].  The  program 
is  based  on  classical  laiminated  plate  theory  and  the  displace- 
ment method  of  finite  element  analysis.  This  progreun  could  be 
used  to  analyze  plies  of  several  different  combinations,  one 
of  which  was  (0,  +45,  90)  , and  included  various  standard  elements, 

5 

including  the  constant  strain  triangle.  The  element  stiffness 
matrix  is  modified  by  either  changing  the  relative  percentage 
of  the  plies  contained  in  the  element  leuninate  or  by  changing 
the  element  material. 


29 


There  are  many  areas  that  should  be  considered  when  con- 
structing a finite  element  model  [32].  These  many  areas  can  be 
condensed  to  four  questions  when  considering  a specific  problem: 


I 

f 

^ i 

i 

I 


I ' 
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1)  What  mesh  model  is  best  for  the  given  problem? 

2)  How  can  symmetry  and  boundary  conditions  be  incorpor- 
ated in  the  model? 

3)  How  accurate  are  the  results  and  has  the  solution  con- 
verged to  the  exact  solution? 

4)  What  do  the  stresses  and  strains  which  result  from  the 
model  actually  mean? 

These  questions  will  now  be  answered  as  they  apply  to  a center 
crack  in  a laminated,  composite  plate. 

In  the  finite  element  analysis  of  cracks  either  standard 
elements  are  used  and  the  mesh  around  the  crack  is  refined  to 
account  for  the  high  stress  gradients,  or  a special  crack  tip 
element  is  used  which  models  the  theoretically  infinite  stress 
at  the  crack  tip.  Since  it  is  not  absolutely  proven  that  all 
ply  stresses  are  infinite  at  the  crack  tip  in  composites,  and 
the  notch  which  simulates  a crack  does  not  extend,  the  standard 
constant  strain  triangle  is  used  in  this  analysis. 

When  elements  such  as  the  constant  strain  triangle  are 
used,  it  is  necessary  to  refine  the  elements  in  the  vicinity 
of  the  crack.  In  one  work  which  studied  the  effect  of  element 
size  [30],  it  was  found  that  accurate  results  could  be  obtain- 
ed by  reducing  the  element  size  in  the  vicinity  of  the  crack 

so  that  the  ratio  of  element  area  to  crack  length  squared 
( \ —6  —6 

I — j was  between  1.2  x 10  and  20  x 10  * These  ratios  re- 

V / 

suited  in  approximately  five  percent  error  in  the  determination 
of  the  isotropic  stress  intensity  factor.  For  this  analysis  two 
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o»K 


finite  element  meshes  were  generated  making  use  of  the  infor- 


mation presented  in  reference  30.  In  the  first/  Figure  13, 

with  163  nodes  and  282  elements,  the  — ratio  was  approximate- 
-6 

ly  20  X 10  . In  the  second.  Figure  14,  with  252  nodes  and  457 

A-.  _g 

elements,  the  — ratio  was  10  x 10  . The  accuracy  of  these 

a* 

two  meshes  when  applied  to  the  composite  crack  problem  was 
assessed  by  checking  stress  convergence. 

Before  these  element  grids  could  be  checked  for  convergence, 
the  symmetry  and  boimdary  conditions  had  to  be  investigated. 
Since  this  is  a two  dimensional  analysis  and  no  bending  loads 
were  applied,  the  Icunina  stacking  sequence  did  not  affect  the 
solution;  therefore,  it  is  possible  to  assume  symmetry  about 
both  the  x-axis  and  the  y-axis.  Due  to  symmetry,  it  is  only 
necessary  to  model  one  quarter  of  the  plate.  The  boundary 
conditions  required  are  to  restrain  the  y-displacements  along 
the  x-axis  and  the  x-displacements  along  the  y-axis  as  shown  in 
Figures  13  and  14. 

The  values  of  were  calculated  for  both  meshes  along  the 
radial  line  running  from  the  crack  tip  parallel  to  the  x-axis 
in  order  to  check  for  convergence.  As  can  be  seen  in  Figure 
15,  the  meshes  yield  values  that  are  essentially  the  same  ex- 
cept in  the  immediate  vicinity  of  the  crack  tip. 

To  check  the  second  mesh  for  accuracy,  the  values  of 
along  the  radial  line  running  from  the  crack  tip  parallel  to 
the  x-axis  were  numerically  integrated,  Fig.  16,  and  the  values 


of  along  the  crack  flank  parallel  to  the  y-axis  were  inte- 
grated, Fig.  17.  From  energy  considerations,  it  is  obvious  that 
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the  Ny  values  should  integrate  to  zero  and  the  values  of 
should  integrate  to  the  applied  load,  P.  The  closeness  of 
these  integration  values  to  the  expected  values  is  a measure 
of  the  accuracy  of  the  numerical  solution.  The  relative 
values  of  N^/P  integrate  to  -.009,  and  the  relative  values  of 
N^/P  integrate  to  .955.  These  results  indicate  that  for  the 
given  assumptions  the  numerical  results  are  slightly  lower  than 
the  actual  stresses  and  that  no  gross  numerical  inaccuracies 
are  introduced. 

A second  convergence  problem  is  the  determination  of  the 
proper  loading  increments.  In  the  progressive  failure  model, 
the  load  displacement  diagraun  becomes  nonlinear.  Fig.  18.  In 
the  nonlinear  portion  of  the  curve,  it  is  necessary  to  reduce 
the  size  of  the  loading  increments  in  order  to  insure  conver- 
gence to  the  correct  solution.  In  this  analysis,  the  load  in- 
crement is  halved  when  the  modelled  displacement  in  a load 
iteration  is  10%  greater  than  the  displacement  detezmined  in 
the  previous  iteration.  The  process  of  reducing  the  increment 
is  continued  until  convergence  is  obtained  over  all  loads  or 
the  slope  of  the  load  displacement  diagram  changes  by  more 
theui  a factor  of  fifty. 

The  last  question  which  needs  to  be  addressed  in  the  finite 

element  modelling  is  how  to  use  the  stresses  which  result  from 

the  analysis.  Since  constant  strain  triangles  are  used,  the 

stress  is  assumed  to  be  constant  over  the  entire  triangle.  For 

the  purpose  of  analyzing  the  stress  fields,  it  is  assumed  that 

the  stress  results  are  for  the  centroid  of  the  triangle.  As 

3: 
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Suggested  in  reference  [32],  values  of  stresses  are  averaged 
over  adjoining  elements  using  the  relative  volume  of  the  el- 
ements as  a weighting  factor.  This  process  reduces  the  varian- 
ce which  exists  where  there  are  high  stress  gradients. 

Bounding  the  Daunage  Zone. 

As  part  of  the  results  from  the  finite  element  analysis, 
the  value  of  the  Tsai-Hill  failure  criterion  for  each  ply  in 
an  element  is  calculated.  Since  the  dcimage  zone  boundary  is 
the  locus  of  points  where  the  value  of  the  failure  criterion 
is  equal  to  one,  it  is  necessary  to  interpolate  between  cal- 
culated values  to  obtain  the  location  of  the  dcunage  zone  bound- 
ary. 

The  interpolation  is  performed  in  the  following  manner. 

The  plate  is  divided  into  fifteen  degree  segments,  where  the 
crack  tip  is  the  origin  of  the  fifteen  degree  radial  lines.  The 
radii  of  the  points  along  these  lines  where  the  values  of  the 
failure  crite.^^ion  are  immediately  above  and  below  one  are 
determined.  Since  the  distances  between  these  points  is  small 
and  the  failure  criterion  varies  in  a monotonic  fashion  along 
the  radials,  the  location  of  the  boundary  is  determined  by  lin- 
early Interpolating  between  the  points  immediately  above  and 
below. 

General  Procedure 

There  are  two  general  procedures  which  are  followed  in  the 


numerical  analysis  stage.  One  is  for  the  purely  elastic  method 

and  the  other  is  for  the  method  which  employs  ply  removal.  In  the 
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elastic  method,  loads  are  applied  in  increments  of  10%  up  to 
the  experimentally  determined  strength.  At  each  level  of  load- 
ing the  bounds  of  the  damage  zone  are  determined,  but  the  ele- 
ment stiffnesses  are  not  changed.  In  the  method  using  ply  re- 
moval, the  seune  general  procedure  is  initially  followed,  but 
after  each  increment  of  loading  the  element  plies  which  failed 
are  removed  before  the  next  increment  of  load  is  applied.  If 
the  displacement  differs  by  more  than  10%  for  the  saune  load  in- 
crement, the  load  increment  is  halved  until  the  model  converges. 
This  procedure  is  continued  until  either  the  experimentally 
determined  strength  is  reached  or  the  model  becomes  unstable, 
where  instability  is  defined  by  the  point  where  the  slope  changes 
by  more  them  a.  factor  of  fifty. 


IV.  Results 


In  this  section,  the  results  of  the  analysis  will  be  presen- 
ted. First,  the  general  load  displacement  and  damage  zone  de- 
termination eure  discussed  followed  by  the  correlation  between 
ply  subcrack  length  and  and  . The  last  section  pertains 
to  fracture  strength  prediction. 

General  Results 

The  load-displacement  data  for  the  elastic  and  progressive 
failure  models  are  given  in  Table  III.  The  load  displacement 
curve  for  the  elastic  model  is  shown  in  Fig.  19.  The  relation- 
ship between  the  load  and  the  displacement  remains  linear. 

The  load  displacement  curves  for  the  progressive  failure 
model  are  shown  in  Fig.  20.  Four  iterations  were  required  be- 
fore the .model  converged.  The  first  iteration  diverged  from 
the  elastic  curve  by  more  than  10%  at  90%  of  the  experimental 
fracture  strength;  this  loading  sequence  was  continued,  though, 
until  an  instability  was  reached  at  110%  of  the  fracture  strength 
The  next  iteration  started  at  85%  of  the  fracture  strength  and 
proceeded  in  5%  increments.  This  iteration  diverged  at  the  90% 
level  again  and  developed  an  instability  at  100%  of  fracture 
strength.  During  the  third  iteration,  divergence  occurred  at 
90%  and  instability  at  92.5%.  The  last  iteration  started  at 
87.5%,  and  developed  an  instability  at  90%. 


Table  III 

Load  Displacement 

Curve  Data 

Per  Cent  of 

Applied 

Applied 

Load 

Experimental 

Load 

Stress 

Displacement 

Fracture  Load 

(lb.) 

(ksi) 

(in. ) 

Elastic  Model 

- 

10 

256 

2.98 

2.1600 

X 

10I3 

20 

504 

5.88 

4.2500 

X 

30 

760 

8.86 

6.4117 

X 

10_^ 

40 

1016 

11.84 

8.4300 

X 

^0-2 

50 

1264 

14.73 

1.0664 

X 

^0-2 

60 

1520 

17.72 

1.2823 

X 

^0-2 

70 

1768 

20.61 

1.4916 

X 

10  , 

80 

2024 

23.59 

1.7075 

X 

^0-2 

90 

2280 

26.58 

1.9235 

X 

10-2 

100 

2528 

29.47 

2.1327 

X 

10 

- 

Progressive 

Failure 

Model, 

1st  Iteration 

-3 

10 

20 

30 

256 

504 

760 

2.98 

5.88 

8.86 

2.1583 

4.2503 

6.4120 

X 

X 

X 

10-3 

10-3 

10-3 

40 

1016 

11.84 

8.5866 

X 

10-2 

50 

1264 

14.73 

1.0705 

X 

10-2 

60 

1520 

17.72 

1.2947 

X 

70 

1768 

20.61 

1.5197 

X 

10.2 

80 

2024 

23.59 

1.7807 

X 

90 

2280 

26.58 

2.1533 

X 

10-2 

100 

2528 

29.47 

3.2329 

X 

10- 

110 

2784 

32.45 

1.2248 

Progressive 

Failure 

Model , 

2nd  Iteration 

-2 

85 

2152 

25.09 

2.0324 

X 

10-2 

90 

2280 

26.58 

2.4415 

X 

10-2 

95 

2408 

28.07 

3.2135 

X 

10  ^ 

100 

2528 

29.47 

1.2352 

Progressive 

Failure 

Model , 

3rd  Iteration 

-2 

87. 

5 

2216 

25.83 

2.3730 

X 

10-2 

90 

2280 

26.58 

3.0427 

X 

10  ^ 

92. 

5 

2344 

27.37 

1.1400 

Progressive 

Failure 

Model , 

4th  Iteration 

-2 

88. 

75 

2248 

26.20 

3.0000 

X 

10 

90 

2280 

26.58 

1.1000 
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Load  Displacement,  5 (in.) 

Figure  19.  Load  Versus  Qisplacement 
Curve  for  Elastic  Model. 
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Load  Displacement,  6 (in.) 

Figure  20.  Load  Versus  Displacement  Curve 
for  Progressive  Failure  "todel. 


The  damage  zone  diagrams  for  the  elastic  model  are  shown 
in  Figures  21-30.  The  plots  correspond  to  the  portion  of  the 
plate  shown  in  the  figures.  The  coordinates  x and  y are 
normalized  by  dividing  them  by  the  plate  width,  w.  The  dotted 
line  parallel  to  the  y-axis  represents  the  crack  and  the  crack 
tip  is  located  at  y/w  =.13.  The  boundary  of  the  plate  is  at 
y/w  =.5.  ' 

Several  generalizations  can  be  drawn  from  examination  of 
these  plots.  First  the  damage  zone  remains  very  small  until 
40%  of  the  fracture  strength  is  reached.  From  Fig.  28,  which 
is  for  80%  of  the  fracture  strength,  the  general  shape  of  the 
damage  zones  in  each  ply  can  be  seen. 

The  damage  zone  in  each  ply  generally  extends  the  farthest 
in  a direction  perpendicular  to  the  fiber  direction.  This  is 
caused  by  the  large  amount  of  shearing  and  transverse  failure. 
The  damage  zone  extends  slightly  behind  the  crack  tip  in  all 
plies. 

The  90“  ply  has  the  largest  damage  zone,  and  the  zone  in 
this  ply  grows  the  fastest.  At  100%  of  the  fracture  load,  all 
of  the  90“  ply  has  failed  except  for  a small  area  along  the 
crack  flanks. 

The  boundaries  for  the  +45“  ply  and  the  -45“  ply  essential.'’.y 
coincide.  The  small  deviances  which  occur  are  believed  to  be 
caused  by  nvimerical  errors.  As  previously  mentioned,  the  stress 
fields  in  the  +45“  plies  showed  the  greatest  disparity  among 
adjoining  elements.  The  damage  zone  in  the  +45“  plies  remains 

essentially^  circular  in  nature. 
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Figure  21.  Damage  Zone  Prediction  at  10%  of  Experimental 
Fracture  Load,  Using  Elastic  Model. 
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Figure  24.  Damage  Zone  Prediction  at  40%  of  EScperimental 
Fracture  Load,  using  Elastic  Model. 
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Figure  27.  Daunage  Zone  Prediction  at  70%  of  Experimental 
Fracture  Load.  Using  Elastic  Model. 
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The  smallest  damage  zone  is  that  for  the  0®  plies.  As  can 


be  seen,  it  is  narrow  and  pointed.  The  area  bounded  by  the  0® 
damage  zone  is  also  the  area  in  which  all  plies  have  failed. 

The  growth  of  this  damage  zone  is  mainly  in  the  direction  of 
the  original  notch.  Due  to  the  shape  and  extension  direction 
of  this  deunage  zone,  the  assumption  that  it  represents  pseudo- 
crack extension  is  validated. 

The  dcimage  zone  boundaries  for  the  progressive  failure  model 
are  shown  in  Figures  31  through  40.  The  boundaries  for  the 
elastic  model  and  the  progressive  failure  model  are  similar  up 
to  the  50%  load  level.  At  this  point,  the  damage  zones  in  the 
progressive  failure  model  begin  to  grow  faster.  Another  dis- 
similarity is  that  the  growth  of  the  damage  zone  behind  the 
crack  tip  stops  at  50%.  The  removal  of  element  plies  causes 
enough  stress  relaxation  so  that  the  stresses  in  this  area  are 
similar  to  the  stresses  along  the  crack  flanks. 

At  the  85%  load  level,  the  +45®  plies  and  the  90®  ply  have 
failed  in  the  entire  region  between  the  crack  tip  and  the  edge 
of  the  plate.  Contrary  to  what  happens  with  the  elastic  model, 
the  90®  ply  never  fails  over  the  entire  plate.  Ultimate  failure 
becomes  imminent  when  0®  ply  damage  zone  reaches  the  edge  of  the 
plate  at  88.75%  of  the  experimental  fracture  strength. 

Several  observations  can  be  made  by  comparing  the  elastic 
model  and  the  progressive  failure  damage  zones.  First,  in  the 
progressive  failure  model,  the  stress  relaxation  caused  by 
element  ply  removal  causes  the  area  where  the  stresses  are  re- 
latively small  to  expand  from  the  immediate  vicinity  of  the 
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Figure  31.  Damage  Zone  Prediction  at  20%  of  Experimental 

Fracture  Load,  using  Progressive  Failure  Model. 
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Figure  32.  Damage  Zone  Prediction  at  30%  of  Experimental 
Fracture  Load,  Using  Progressive  Failure  Model. 
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. Figure  34.  Damage  Zone  Prediction  at  50%  of  Experimental 

j Fracture  Load,  Using  Progressive  Failure  Model. 
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Figure  37.  Damage  zone  Prediction  at  80%  of  Experimental 
Fracture  Load,  Using  Progressive  Failure  Model 
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Figure  39.  namage  Zone  Prediction  at  87.5%  of  Experimental 
Fracture  Load,  Using  Progressive  Failure  Model. 
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crack  flanks,  which  is  where  the  relatively  small  stresses 
occur  in  the  elastic  model.  Next,  as  is  obvious  by  com- 
paring respective  damage  zone  diagrams,  the  damage  zones 
predicted  using  the  progressive  failure  model  are  larger 
than  those  predicted  using  the  elastic  model.  As  a conse- 
quence of  this  larger  size,  the  predicted  failure  load  is 
less  for  the  progressive  failure  model  than  it  is  for  the 
elastic  model.  Last,  it  must  be  recognized  that  these  two 
models  are  intended  to  boiind  the  actual  case.  The  actual 
d2unage  zone  boundaries  will  exist  somewhere  between  those 
predicted  with  the  progressive  failure  model  and  those 
predicted  with  the  elastic  model. 

Correlation  of  Subcrack  Length 

The  first  application  of  -the  finite  element  analysis  is 
to  check  the  correlation  between  the  subcrack  length  in  each 
ply  and  the  values  of  and>C  . Two  values  of  are  cal- 
culated; one  is  for  the  purely  elastic  analysis,  and  the 
other  is  for  the  progressive  failure  analysis. 

The  first  task  is  to  determine  the  length  of  the  subcracks 

in  each  ply.  The  subcrack  measurements  are  obtained  from  the 

photographs  shown  in  Figures  4 through  9.  As  can  be  seen  in 

these  photographs  the  exact  length  of  the  ply  subcracks  is  not 

easily  discernible.  The  subcracks  perpendicular  to  the  notch 

are  in  the  0®  ply;  the  subcracks  colinear  with  the  notch  are 

in  the  90®  ply,  and  the  subcracks  running  oblique  to  the  notch 

are  in  the  +45®  plies.  Since  the  orientation  of  the  plate  is 
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unknown,  it  is  impossible  to  discern  between  the  +45**  and 
-45*  plies.  The  subcrack  lengths  shown  in  Table  IV  are  the 
average  of  the  measured  lengths. 


Table  IV 

Subcrack  Lengths 

Applied 

0* 

90** 

+45** 

Stress 

ply 

ply 

ply 

(ksi) 

(in.) 

(in.) 

(in.) 

14.734 

.031 

.047 

.023 

17.718 

.055 

.086 

.039 

20.609 

.057 

.115 

.046 

23.593 

.060 

.226 

.068 

26.577 

.084 

.436 

.092 

28.069 

.087 

.564 

.095 

The  values  of  are  obtained  from  equation  (4);  and 
is  calculated  using  equation  (8) . The  crack  half  length, 
a,  is  calculated  using  equation  (5) . The  computed  values  of 
Kj  at  the  various  stress  levels  for  the  elastic  analysis  are 
shown  in  Table  V. 

The  values  of  for  the  progressive  failure  analysis 
are  shown  in  Table  VI.  This  table  does  not  go  to  the  same 
stress  level  as  the  elastic  analysis  since  the  analysis  de- 
veloped an  instability  before  100%  of  the  experimental 
notched  strength  was  reached.  The  values  of  ^ are  shown  in 

Tcdsle  VII.  The  compliance  values  are  obtained  from  Fig.  20. 

2 

The  values  of  determined  from  the  elastic  analysis 

versus  the  subcrack  length  in  each  ply  are  shown  in  Fig.  41. 

As  was  found  in  reference  4,  there  anoears  to  be  a linear 
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Table  V 

Calculations 

(Elastic  Model) 

Applied 

Crack 

Plate 

*^1 

Stress 

Length 

Width 

L. 

(ksi) 

(in. ) 

(in.) 

(ksi-in. ) 

(ksi-in.) 

2.9841 

.2559 

1.9685 

2.75 

7.59 

5.8750 

.2559 

1.9685 

5.00 

25.04 

8.8591 

.2640 

1.9685 

8.32 

69.24 

11.8433 

.2687 

1.9685 

11.04 

121.85 

14.7341 

.2745 

1.9685 

14.15 

200.19 

17.7183 

.2814 

1.9685 

17.28 

298.44 

20.6091 

.2865 

1.9685 

20.28 

411.36 

23.5933 

.2928 

1.9685 

23.51 

552.87 

26.5774 

.3034 

1.9685 

27.04 

731.34 

29.4683 

.3116 

1.9685 

30.46 

927.85 

Table  VI 

Calculations 

(Progressive  Failure  Model) 

Applied 

Crack 

Plate 

■'i' 

Stress 

Length 

VJidth 

u 

9 

(ksi) 

(in. ) 

(in. ) 

(ksi-in. ) 

(ksi-in. ) 

2.9841 

. 2559 

1.9685 

2.75 

7.59 

5.  8750 

.2559 

1.9685 

5.00 

25.04 

8.8591  ' 

.2615 

1.9685 

8.28 

68.50 

11.8433 

.2701 

1.9685 

11.27 

126.98 

14.7341 

.2881 

1.9685 

14.55 

211.62 

17.7183 

.3213 

1.9685 

18.65 

344.93 

20.6091 

.3769 

1.9685 

23.96 

573.85 

23.5933 

.4487 

1.9605 

30.89 

953.51 

25.0853 

.5698 

1.9685 

39.89 

1591.36 

25.8341 

.4627 

1.9685 

59.66 

3558.77  f- 
16.9  X 10 

26.2044 

« 

.9842 

1.9685 

4115.75 
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Table  VII 

strain  Energy  Release  Rates 

Squared 

Change  in 

Crack 

Strain  Energy 

Lo^d 

Compliance 

Extension 

Release  Rate 

P* 

Ac 

Aa 

/ ^ \ 

(Ib^) 

(in. /lb.) 

(in. ) 

/ in.-lb.l 

\ in.  / 

1597696 

1.77  x 10”® 
4.87  X 10"“ 

7.78  X 10”? 

3.60  X lO”^ 

.39 

2310400 

6.64  X 10”i 

.85 

3125824 

1.11  X 10"t 

1.09 

4096576 

2.02  X lO”, 

1.44  X 10"7 

2.89 

4631104 

6.46  X lO”' 

2.42  X 10"7 

6.18 

4910656 

5505354 

1.26  X 10”° 
2.64  X 10”° 

3.86  X 10“r 
4.43  X lO"-^ 

8.05 

15.04 

relation  between  and  the  subcrack  lengths. 

2 

The  values  of  determined  from  the  progressive  failure 

analysis  versus  subcrack  lengths  are  shown  in  Fig.  42.  For  the 
2 

values  of  , considering  values  of  load  from  14.7341  ksi  to 

23.593  ksi,  the  correlation  appears  to  be  linear.  At  the  last 

load  level  where  the  damage  zone  extends  to  the  edge  of  the  plate, 

there  is  no  correlation  between  and  the  subcrack  lengths. 

The  values  of  ^ versus  subcrack  length  are  plotted  in  Fig. 

2 

43.  As  with  the  versus  subcrack  length  for  the  progressive 
failure  case,  the  relation  appears  to  be  linear.  The  data  point 
associated  with  the  last  load  increment  is  not  within  the  range 
of  a linear  relation  since  the  equation  used  to  approximate  the 
derivative  of  structural  compliance  with  respect  to  crack  exten- 
sion, equation  (7)  , is  not  valid. 

From  the  limited  amount  of  experimental  data,  it  is  difficult 

^•o  determine  absolutely  if  a linear  relation  exists.  For  the 

2 M 

data  available  though,  values  of  and  A-  are  linearly  related 


flqur*  41.  Versus  Subcrack  Length 

for  Elastic  Model. 


2 

Figure  42.  K,  Versus  Subcrack  Length 

for  Progressive  failure  Model. 
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to  the  ply  subcrack  lengths  except  when  the  load  is  at  the 
point  where  the  model  exhibits  large  nonlinear  behavior. 

Failure  Prediction 

As  was  stated  in  the  theory  chapter,  calculated  values  of 
the  opening  mode  stress  intensity  factor,  K^,  can  be  used  to 
predict  failure.  For  this  specimen  the  critical  stress  in- 
tensity  factor  is  between  40.22  and  41.67  ksi-in.  Referring 
to  Table  V,  it  is  seen  that  failure  would  occur  at  some  value 
over  29.4683  ksi  which  is  the  experimental  notched  strength. 

If  the  value  of  continued  to  increase  in  the  same  manner, 
the  predicted  fracture  strength  would  be  approximately  39  ksi 
or  32%  over  the  experimental  strength.  Using  the  values  of 
Kj  in  Table  VI,  the  failure  strength  would  be  between  25.0853 
and  25.8341  ksi.  This  is  in  error  by  12-15%.  As  was  expected 
the  elastic  analysis  provides  an  upper  bound  on  the  fracture 
strength  and  the  progressive  failure  analysis  provides  a lower 
bound . 

The  load-displacement  diagram.  Fig.  20,  can  be  used  to 
determine  the  fracture  load  from  a stability  standpoint.  After 
the  last  iteration  the  slope  changes  from  5.104  x 10^  lb. /in., 
for  the  load  increment  from  87.5%  to  88.75%  of  the  experimental 
fracture  load,  to  2.99  x 10^  lb. /in.  for  the  increment  from 
88.75%  to  90%.  Therefore,  the  failure  strength  becomes  the 
stress  at  88.75%  or  26.204  ksi.  This  is  below  the  actual  frac- 
ture strength  by  11%.  For  the  instability  analysis,  there  is 

not  an  upper  bound  since  the  elastic  load  displacement  diagram 
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remains  linear. 

The  last  method  that  can  be  used  to  predict  fracture 
strength  is  the  applied  load  versus  load  bearing  area  diagram 
(P-Aj^g  diagram) . The  values  of  the  applied  load  and  the  re- 
maining load  bearing  area  data  are  shown  in  Table  VIII. 


Table  VIII 


Load  and  Load  Bearing  Area  Data 


Load,P 

(lb.) 


256 

504 

760 

1016 

1264 

1520 

1768 

2024 

2152 

2216 

2248 

2280 

2528 


Elastic  Remaining 
Load  Bearing  Area 
Alb (in.*) 


6.35  X 10"^ 
6.35  X lO", 
6.28  X lO", 
6.24  X lO", 
6.19  X 10"; 
6.12  X 10", 
6.08  X 10"; 
6.03  X 10"“^ 


5.93  X 10 
5.86  X 10' 


Progressive  Failure 
Remaining  Load 
Bearing  Area  A^g(in, 


6.35 

6.35 

6.30 

6.22 

6.07 

5.78 

5.29 

4.67 

3.61 

1.93 


X 10 
X 10 
X 10 
X 10 
X 10 
X 10 
X 10 
X 10 
X 10 
X lO' 
0 


* Analysis  was  not  performed  for  these  loads. 


The  P versus  A^^g  diagram  is  shown  in  Fig.  44.  The  values 

along  the  horizontal  axis  correspond  to  values  of  A^g  which  is 

the  load  bearing  area  of  the  plate  between  the  notch  and  the 

edge  of  the  plate.  The  load  bearing  area  represents  the 

portion  of  the  plate  in  which  all  lamina  have  not  failed.  The 

vertical  axis  values  are  the  applied  loads,  P.  The  straight 

line  running  in  an  oblique  direction  from  the  origin  represents 

the  boundary  between  loads  and  load  bearing  areas  which  do  not 
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The  slope  of  this  line  is  equal  to  the  unnotched  failure 
stress  of  54.4  ksi. 

In  order  to  explain  the  significance  of  the  points  on  this  1 

i 

diagram,  the  effect  of  notch  sensitivity  is  ex^unined.  If  the 
plate  is  not  notch  sensitive,  the  plate  would  fail  when  the 
load  per  area  exceeded  the  notched  tensile  stress.  The  failure 
line  would  extend  parallel  to  the  load  axis  and  the  failure 
load  would  be  3400  lb.  Since  the  plate  is  notch  sensitive,  the 
failure  load  is  less,  2530  lb.  This  is  a 34%  error. 

The  growth  of  a damage  zone  at  the  crack  tip  accounts  for 
the  notch  sensitivity.  The  elastic  analysis  can  be  used  to 
model  damage  zone  growth  as  shown  in  the  diagram.  Using  only 
the  elastic  analysis,  the  predicted  fracture  load  is  3100  lb. 

This  is  a 23%  error.  As  expected,  the  elastic  analysis  pro- 
vides a prediction  which  is  above  the  actual  fracture  load. 

The  progressive  failure  analysis  can  also  be  used  to  predict 
fracture  strength.  The  failure  curve  using  this  analysis  be^- 
comes  nonlinear  in  the  upper  load  levels  as  the  damage  zone 
growth  accelerates.  The  predicted  fracture  load  using  the  pro- 
gressive failure  analysis  is  2110  lb.  This  prediction  is  16% 
below  the  actual  fracture  strength.  As  predicted  using  only 
theoretical  considerations,  the  progressive  failure  analysis 
provides  a lower  bound  on  the  fracture  strength. 


] 


77 


V.  Conclusions 

As  cem  be  seen  from  the  dcunage  zone  diagreuns,  the  eunount 
of  damage  in  each  ply  is  vastly  different.  Since  it  is  impos- 
sible to  experimentally  measure  this  damage  zone,  the  use  of 
a numerical  model  such  as  the  one  presented  in  this  thesis  is 
warranted.  Of  course  this  thesis  only  studied  one  launinate 
with  one  notch  orientation,  but  the  accuracy  of  the  model 
would  indicate  that  further  study  should  be  conducted. 

Using  numerical  models,  it  was  possible  to  correlate  ply 
subcrack  length  and  two  fracture  mechanics  pareuneter.  Although 
it  was  not  shown  that  the  subcrack  length  was  related  to 
ultimate  failure,  the  linear  relation  between  , and 

subcrack  length  does  indicate  that  some  principles  of  fracture 
mechanics  do  apply,  at  least  in  the  immediate  vicinity  of  the 
crack  tip'. 

Through  numerical  modelling,  it  was  also  possible  to  bound 
the  fracture  strength  using  either  or  one  of  the  other  models. 
Although  the  instability  analysis  provided  the  closest  approxi- 
mation of  fracture  strength,  it  did  not  provide  an  uppe..  bound. 
The  approach  provided  bounds  which  were  as  close  to  the 

actual  strength  as  that  calculated  using  values.  Since  the 
P-Aj^g  approach  could  be  applied  to  all  types  of  laminates  and 
structures,  it  is  considered  better. 

The  numerical  model  presented  in  this  paper  is  crude  and 

could  obviously  be  improved  uoon.  The  most  important  area 
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requiring  improvement  is  that  of  the  strength  criterion.  A 
better  criterion  could  possibly  give  better  estimates  of  which 
element  plies  have  failed.  The  next  area  requiring  improve- 
ment is  in  the  treatment  of  element  stiffness  after  failure. 
After  exceeding  the  failure  criterion,  the  element  probably 
retains  some  load  carrying  capability.  Since  this  analysis 
completely  discounted  all  stiffness  after  element  ply  failure, 
it  should  provide  estimates  which  are  conservative.  Yet  the 
method  always  provides  a lower  bound  solution  which  is  impor- 
tant when  considering  problems  in  which  experimental  data  is 
nonexistent  or  the  experiment  is  in  the  planning  stage.  Im- 
proving stiffness  characteristics  may  not  provide  a closer, 
conservative  result.  The  last  procedure  which  requires  im- 
provement is  in  the  method  of  loading.  Through  the  use  of 
more  sophisticated  incremental  loading  and  convergence  methods, 
associated  with  nonlinear  analysis,  the  predictions  could  be 
improved . 

The  applicability  of  the  finite  element  method  in  analyzing 
composite  fracture  has  been  shown  for  this  special  case.  Since 
finite  element  models  can  be  applied  to  complicated  structures, 
and  are  not  as  costly  as  experimentation,  further  study  into 
the  application  of  finite  elements  to  this  type  of  problem  is 
necessary  and  should  prove  profitable. 
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